### Analysis of results for simulations for the following case: admit pivot = 0.67; decision pivot = 0.5
### Generates Figure 1b and Table A2 in manuscript
### Last Revised and commented: 9 June 2016, JBS

##########################################################################################

# clear workspace
rm(list=ls())

# load libraries
library(reldist)
library(memisc)
library(Hmisc)

# set path (needs to be changed by user); path will depend on where the folder is stored within in the file structure
path <- "~/Dropbox/IO imperfect monitoring/GLS_II_REP"

# set working directory based on path. Should work automatically if path is specified correctly
setwd(path)

# load initial membership data
load('initmems.Rdata')

# set working director to folder, where the results are stored
setwd(paste(path , "a67d50" , sep = "/"))

# set number of simulations equal to the number of stored simulations
sims<-300

# create a holding bin for results from the simulations
memsims <- vector("list",length=sims)

# read in results from simulations
for (i in 1:sims){
	
	inputFile <- paste("res_",i,".csv",sep="")
	con  <- file(inputFile, open = "r")
	
	dtalist <- vector("list", length=100)
	
	counter <- 1
	
	while (length(oneLine <- readLines(con, n = 1, warn = FALSE)) > 0) {
    	 
    	 dtalist[[counter]] <- as.numeric(unlist(strsplit(oneLine, " ")))
         counter <- counter+1
         } 
	close(con)
	memsims[[i]]<-dtalist
}

# calculate median of initial memberships
initmed <- unlist(lapply(samp,median))

# calculate range of initial memberships
initrange <- unlist(lapply(samp,max))-unlist(lapply(samp,min))

# calculate Gini of initial memberships
initgini <- unlist(lapply(samp,gini))

# calculate size of initial memberships
initsize <- unlist(lapply(samp,length))

# create holding bins for changes in median/Gini/size between initial membership and final simulated membership
medianP <- vector("list", length= sims)
medianG <- vector("list", length= sims)
medianS <- vector("list", length= sims)

# calculate and store changes in median/Gini/size between initial membership and final simulated membership
for (i in 1: sims) {
	
	medianP[[i]] <- abs(initmed - unlist(lapply(memsims[[i]], median)))
	
	medianG[[i]] <-  unlist(lapply(memsims[[i]], gini)) / initgini 
	
	medianS[[i]] <- unlist(lapply(memsims[[i]], length)) - initsize 

}

# putting change calculations in data frames
medianP <- as.data.frame(do.call(cbind,lapply(medianP,unlist)))
mmedP <- apply(medianP,1,median) 
medianG <- as.data.frame(do.call(cbind,lapply(medianG,unlist)))
mmedG <- apply(medianG,1,median) 
medianS <- as.data.frame(do.call(cbind,lapply(medianS,unlist)))
mmedS <- apply(medianS,1,median) 

# Standardized regressions to replicate Table A2 in the manuscript

r1a <- lm(scale(mmedP) ~ scale(initsize)+scale(initrange)+scale(initgini))
cf.r1a <- coef(summary(r1a))[-1,1]
se.r1a <- coef(summary(r1a))[-1,2] 

r2a <- lm(scale(mmedG) ~ scale(initsize)+scale(initrange)+scale(initgini))
cf.r2a <- coef(summary(r2a))[-1,1]
se.r2a <- coef(summary(r2a))[-1,2] 

r3a <- lm(scale(mmedS) ~ scale(initsize)+scale(initrange)+scale(initgini))
cf.r3a <- coef(summary(r3a))[-1,1]
se.r3a <- coef(summary(r3a))[-1,2] 

# LaTeX code for standardized regressions to replicate Table A2 in the manuscript
toLatex(mtable(r1a,r2a,r3a))

#### Code for replicating Figure 1(b) in manuscript 

var.names <- rep(c("Initial Size","Initial Range","Initial Dispersion"),3)

yaxis1 <- c(3.1,3,2.90)
yaxis2 <- c(2.1,2,1.9)
yaxis3 <- c(1.1,1,0.90)

# Put figure in main directory
setwd(path)

pdf("coef_figa67d50.pdf", height = 8, width = 5.25)
#layout(matrix(c(2,1),1,2), #in order to add variable categories and braces to left side of plot, 
#   widths = c(1.7, 5))#we use layout command, create a small second panel on left side.
                #using c(2,1) in matrix command tells R to create right panel 1st
			#layout.show(2) #can use this command to check results of layout command (but it must be commented out when creating PDF).


par(mar=c(2,8,.5,1), lheight = .8)#set margins for regression plot
plot(cf.r1a, yaxis1, type = "p", axes = F, xlab = "", ylab = "", pch = 15, cex = 1.2, #plot model 1 coefs using black points (pch = 19, default = black), adding the "adjust amount" to the y.axis indicator to move points up
    xlim = c(min((cf.r1a-qnorm(.975)*se.r1a -.1), (cf.r2a-qnorm(.975)*se.r2a -.1), (cf.r3a-qnorm(.975)*se.r3a -.1), na.rm = T), #set xlims at mins and maximums (from both models) of confidence intervals, plus .1 to leave room at ends of plots
    max((cf.r1a+qnorm(.975)*se.r1a +.1), (cf.r2a+qnorm(.975)*se.r2a +.1), (cf.r3a+qnorm(.975)*se.r3a +.1), na.rm = T)),  #use na.rm=T since vectors have missing values
    ylim = c(.75,3.25),
     main = "")
axis(1,at = seq(-1,1, by = 0.5), label = seq(-1,1, by = 0.5), mgp = c(3,1,0), cex.axis = .8)#add x-axis and labels; "pretty" creates a sequence of  equally spaced nice values that cover the range of the values in 'x'-- in this case, integers
axis(2, at = c(yaxis1,yaxis2,yaxis3), label = var.names, las = 1, tick = T, cex.axis =1.1)#add y-axis and labels; las = 1 makes labels perpendicular to y-axis


#axis(3,pretty(coef.vec.1, 3))#same as x-axis, but on top axis 
abline(h = c(yaxis1,yaxis2,yaxis3), lty = 2, lwd = .5, col = "light grey")#draw light dotted line at each variable for dotplot effect
#box(bty="l")#draw box around plot
segments(cf.r1a-qnorm(.975)*se.r1a, yaxis1, cf.r1a+qnorm(.975)*se.r1a, yaxis1, lwd =  1.3)#draw lines connecting 95% confidence intervals
#for tick marks indicating 90% ci's use following 2 lines:
#segments(coef.vec.1-qnorm(.95)*se.vec.1, y.axis+adjust -.035, coef.vec.1-qnorm(.95)*se.vec.1, y.axis+adjust +.035, lwd = 1.1)#draw vertical tick marks at 90% 
#segments(coef.vec.1+qnorm(.95)*se.vec.1, y.axis+adjust -.035, coef.vec.1+qnorm(.95)*se.vec.1, y.axis+adjust +.035, lwd = 1.1)   #confidence intervals
abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing

#add 2nd model
    #because we are using white points and do want the lines to go "through" points rather than over them
        #we draw lines first and the overlay points
segments(cf.r2a-qnorm(.975)*se.r2a, yaxis2, cf.r2a+qnorm(.975)*se.r2a, yaxis2, lwd =  1.3)#draw lines connecting 95% confidence intervals
#for tick marks indicating 90% ci's use following 2 lines:
#segments(coef.vec.2-qnorm(.95)*se.vec.2, y.axis-adjust -.035, coef.vec.2-qnorm(.95)*se.vec.2, y.axis-adjust +.035, lwd = 1.1)#draw vertical tick marks at 90% 
#segments(coef.vec.2+qnorm(.95)*se.vec.2, y.axis-adjust -.035, coef.vec.2+qnorm(.95)*se.vec.2, y.axis-adjust +.035, lwd = 1.1)  #confidence intervals
points(cf.r2a, yaxis2, pch = 16, cex = 1.2) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color

#add 3rd model
    #because we are using white points and do want the lines to go "through" points rather than over them
        #we draw lines first and the overlay points
segments(cf.r3a-qnorm(.975)*se.r3a, yaxis3, cf.r3a+qnorm(.975)*se.r3a, yaxis3, lwd =  1.3)#draw lines connecting 95% confidence intervals
#for tick marks indicating 90% ci's use following 2 lines:
#segments(coef.vec.2-qnorm(.95)*se.vec.2, y.axis-adjust -.035, coef.vec.2-qnorm(.95)*se.vec.2, y.axis-adjust +.035, lwd = 1.1)#draw vertical tick marks at 90% 
#segments(coef.vec.2+qnorm(.95)*se.vec.2, y.axis-adjust -.035, coef.vec.2+qnorm(.95)*se.vec.2, y.axis-adjust +.035, lwd = 1.1)  #confidence intervals
points(cf.r3a, yaxis3, pch = 17, cex = 1.2 ) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color

#Create Variable Categories and Braces to go in 2nd plot

#par(mar=c(2,0,.5,0)) #set margins--- bottom (1st number) and top (3rd number) must be the same as in 1st plot
#plot(seq(0,1,length=length(var.names)), c(yaxis1,yaxis2,yaxis3), type = "n", axes = F,  xlab = "", ylab = "")#call empty plot using type="n"
    #use a sequence of length 20 so that x and y have same length

#text(.55, 2.87, "Median Change", font = 3, cex  = 1.1)#Add text; "srt" rotates to 90 degrees, font = 3 == italics
#text(.55, 2, "Dispersion Change", font = 3, cex  = 1.1)#Add text; "srt" rotates to 90 degrees, font = 3 == italics
#text(.55, 1.12, "Expansion", font = 3, cex  = 1.1)#Add text; "srt" rotates to 90 degrees, font = 3 == italics
dev.off()